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The stability of two-dimensional (2D) layers and membranes is subject of a 
long standing theoretical debate. According to the so called Mermin- Wagner 
theorem [I], long wavelength fluctuations destroy the long-range order for 2D 
crystals. Similarly, 2D membranes embedded in a 3D space have a tendency 
to be crumpled [2J. These dangerous fluctuations can, however, be suppressed 
by anharmonic coupling between bending and stretching modes making that a 
two-dimensional membrane can exist but should present strong height fluctua- 
tions [H IHl H]- The discovery of graphene, the first truly 2D crystal [5J, E] and the 
recent experimental observation of ripples in freely hanging graphene |7J makes 
these issues especially important. Beside the academic interest, understanding 
the mechanisms of stability of graphene is crucial for understanding electronic 
transport in this material that is attracting so much interest for its unusual 
Dirac spectrum and electronic properties [81 [91 I10|, 111] , Here we address the 
nature of these height fluctuations by means of straightforward atomistic Monte 
Carlo simulations based on a very accurate many-body interatomic potential for 
carbon |12| . We find that ripples spontaneously appear due to thermal fluctu- 
ations with a size distribution peaked around 70A which is compatible with 
experimental findings [7] (50-100 A) but not with the current understanding of 
stability of flexible membranes [21|3l|4]. This unexpected result seems to be due 
to the multiplicity of chemical bonding in carbon. 

The phenomenological theories for flexible membranes [21 EH E] are derived in the contin- 
uum limit without including any microscopic feature and their applicability to graphene in 
the interesting range of temperatures, sample sizes etc. is not evident. We present Monte 
Carlo simulations of the equilibrium structure of single layer graphene. By monitoring the 
normal-normal correlation functions we can directly compare our results to the predictions of 
the existing theories. The effect of interest is crucially dependent on acoustical phonons and 
interactions between them and therefore the simulations require samples much larger than 
interatomic distances, typically many thousands of atoms at thermal equilibrium, making 
prohibitive ab-initio simulations a la Car-Parrinello [13]. However, for carbon a very accurate 
description of energetic and thermodynamic properties of different phases is provided by the 
effective many body potential LCBOPII [T21 E] . This bond order potential is constructed 
in such a way as to provide a unified description of the energetics and elastic constants of all 
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carbon phases as well as the energy characteristics of different defects with accuracy com- 
parable to experimental accuracy. We find clear deviations from harmonic behavior for long 
wavelength fluctuations but, instead of the expected power-law scaling at long wavelengths, 
we find a marked maximum of fluctuations with wavelength of about 70 A. Against the 
expectations, we also find a stiffening of the bending rigidity with increasing temperature. 
We relate these features to fluctuations in bond length that in carbon signal a partial change 
from conjugated to single/double bonds, with consequent deviations from planarity. 

We perform atomistic Monte Carlo (MC) simulations based on the standard Metropolis 
algorithm for approximately squared samples of different sizes (see Table I) with periodic 
boundary conditions. We always start with completely flat graphene layers to avoid any bias. 
Typically we used 400000 MC steps (1 MC step corresponds to N attempts to a coordinate 
change) to equilibrate and half a million for averaging. Every 5 MC steps, we allow for 
isotropic size fluctuations. The energy of a given configuration is evaluated according to the 
bond order potential LCBOPII developed in the recent past [TJ]. This potential is based 
on a large database of experimental and theoretical data for molecules and solids and has 
been proven to describe very well thermodynamic and structural properties of all phases 
of carbon and its phase diagram in a wide range of temperatures and pressures [121 IH]- 
We believe that LCBOPII can give a good description of graphene because it reproduces 
correctly the elastic properties of graphene and yields structure and energetics of vacancies 
in graphene in good agreement with ab initio calculations [JS]. Of particular importance 
here, is that bond order potentials correlate coordination, bond length and bond strength, 
allowing changes between single, double and conjugated bonds with the correct energetics. 
Most simulations have been performed at room temperature T = 300 K but we have also 
simulated some of the structures at very high temperature T = 3500 K close to the bulk 
graphite melting to study qualitatively temperature effects. 

A typical snapshot of graphene at room temperature is shown in Fig{T] The first thing 
to notice is that height fluctuations are present at equilibrium. We find a very broad dis- 
tribution of height displacements h with a typical size of order of 0.6 A for the N = 8640 
sample, comparable to the interatomic distance 1.42 A. 

The natural way to analyse further our results is to compare them to the results and 
predictions of the phenomenological theories of thermal fluctuations in flexible membranes (21 
[31 0J. To this purpose it is worth to review in some detail the main ideas and results of 
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these theories [2j [3j H] that are not of common knowledge outside the soft condensed matter 
community. The primary quantities are the two-component displacement vector in the plane 
u, the out of plane displacement out of plane h and the normal unit vector n with in-plane 
component — V/i/ yjl + (V/i) 2 illustrated in Figj^J The elastic energy of the membrane is 
given by 

(1) 



E = J d 2 x ^(V 2 h) 2 + fiu 2 af3 + -ul 



where k is the bending rigidity, // and A are Lame parameters and u a p is the deformation 
tensor: 

_ If du a dup dh dh \ 
a/3 2 \dxfs dx a dxadxpj' 
In harmonic approximation, by neglecting the last, non linear, term in the deformation 

tensor (|2]), the bending (h) and stretching (u) modes are decoupled. In this approximation 

the Fourier components of the bending correlation function with wavevector q is 

l\h l 2 \ TN (i\ 

where N is the number of atoms, = L x L y /N is the area per atom and T is the temperature 
in units of energy. In this approximation, the mean square displacement in the direction 
normal to the layer is 

W = £<i^iV^ 2 (4) 

q 

where L is a typical linear sample size. The correlation function of the normals 

G(g) = <|n q | 2 )=g 2 <|/; q | 2 > (5) 

in this approximation becomes 

TN 

Go (q) = ~^— 2 (6) 
KS q 2 

which implies that the mean square angle between the normals is logarithmically divergent 
as L — ► oo [2]. This behaviour indicates the tendency to crumpling of membranes due to 
thermal fluctuations. 

Deviations from this harmonic behaviour, namely anharmonic coupling between bending 
and stretching modes, can stabilize the flat phase by suppressing the long wavelength fluc- 
tuations [21 El H] . In this case the corresponding correlation function of the normals is given 
by the Dyson equation 

G- 1 (q) = G 1 (q) + S (g) (7) 



with self energy 

where go — 2tt^B/k, B being the two-dimensional bulk modulus, rj is the anomalous rigidity 
exponent and A is a numerical factor. 

The simplest way to derive this expression is to use the self-consistent perturbation the- 
ory [4 J which gives 77 ~ 0.8 in reasonable agreement with the results of Monte Carlo sim- 
ulation for a model of tethered membranes [16]. However, Eqs.(|7|,([8j can be written from 
a more general scaling consideration [3], q being the only factor with dimension of an in- 
verse length that can be constructed from relevant parameters of this theory. As a result 
of this anharmonic coupling, the typical height of fluctuations in the direction normal to 
the membrane is much smaller than the one given by Eq.Q and scales with the sample 
size as L*», with ( = 1 — rj/2. Nevertheless, the fluctuations are still anomalously large and 
they can be much larger than the interatomic distance for large samples. Thus, the theory 
predicts an intrinsic tendency to ripple formation. At the same time, the amplitude h oc Lf- 
of these transverse fluctuations remains much smaller than the sample size and preserves the 
long-range order of the normals so that the membrane can be considered as approximately 
flat and not crumpled. 

Another structural issue is the existence of long range crystallographic order in mem- 
branes that can be destroyed by a finite concentration of topological defects, namely dis- 
locations and disclinations. For the case of 2D crystals in 2D space both types of defects 
have infinite energy as L —>■ 00: the elastic energy of dislocations grows as ln(L), and of 
disclinations as L. A general analysis [2] (Chap. 6) for flexible membranes shows that 
these divergencies are suppressed by bending in such a way that the energy of disclina- 
tions behaves as ln(L) whereas that of dislocations remains finite, and of the order of n. 
This means that the orientational order survives whereas translational order is destroyed 
by spontaneous creation of dislocations. However, the corresponding correlation length is 
oc exp(const ■ k/T) which makes this mechanism completely irrelevant for covalently bonded 
layers such as graphene, with k »s 1.17 eV [17] . Indeed, we never observe any topological 
defect in our simulations also at very high temperature, nor any experimental evidence of 
their existence has been reported [7J. 

To obtain a quantitative comparison with these theoretical predictions for the spatial 
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distribution of the ripples, we have calculated numerically the Fourier components of the 
correlation function of the normals G(q) for q x and q y multiples of 2tc/L x and 27T /L y respec- 
tively. To our surprise, our numerical results are not described by this general theory. 

In Fig|3]we show G(q) at T = 300 K for all considered samples and in Fig|4]we compare 
the results at two temperatures for the sample N = 8640. First of all, there is a whole range 
of wave vectors where the harmonic approximation given by Eq.Q is quite accurate. The 
interval is restricted above close to the Bragg peaks and in the limit of small q. The deviations 
are opposite at the two sides. The rigidity k = 1.1 eV extracted from the data at T = 300 K 
by comparison with Eq.(|6]) is in very good agreement with the experimental value of 1.2 eV 
derived from the phonon spectrum of graphite [T7j. Surprisingly, the same comparison at 
T = 3500 K gives a higher value k rs 2.0 eV whereas the contimuum theory [2] predicts the 
opposite trend k(T) ~ k — (3T/47r) In (L/a). This is due to the peculiar character of bonding 
in carbon. In the ground state of graphene all bonds are equivalent. However, even at room 
temperature, there is a large probability of having an asymmetric distribution of short/long 
(strong/weak) bonds, associated with local deviations from planarity. Indeed, the radial 
distribution functions shown in Fig|5] for both temperatures show a broad distribution of 
first neighbors bond lengths, going down to the length of double bonds in carbon of 1.31 
A at high temperature. Changes of bond conjugation are also the reason for the negative 
thermal expansion coefficient in graphene found in ab-initio calculations [18]. We observe a 
0.13 % contraction of the lattice spacing at T = 300 K, in good agrement with the 0.11 % 
of Ref. [TSJ We believe that it is the ability of carbon to form different types of bonding that 
makes graphene different from a generic two-dimensional crystal. 

At small q, the behaviour of G(q) is not described by the harmonic approximation Go(q) 
nor by the anharmonic expression G a (q). The most remarkable feature of G(q) is a maximum 
instead of the power law dependence G a (q) that implies the absence of any relevant length 
scale in the system. The presence of this maximum, instead, means that there is a preferred 
average value of about 70 A. This length is also recognizable in real space images, as shown 
by the arrows in Fig|T} Indeed, the two samples that are smaller than this length do not 
show this decrease of G(q) at low q. The results at T = 3500 K confirm this picture but 
with the maximum shifted to a larger q corresponding to a length of roughly 30 A. This 
temperature dependence of the typical ripple lengthscale should be measurable. 

The preferable length that we find is reminiscent of the "avoided criticality" scenario in 
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frustrated systems near second-order phase transitions [19]. Such systems have an intrinsic 
tendency to be modulated and to form some inhomogeneous patterns that destroys the 
scaling. It is known that some soft condensed matter membranes tend to spontaneous 
bending and ripple formation [201 EU E2] but this behaviour for elemental solids like graphene 
is rather unexpected. 

The results obtained are relevant not only for a better understanding of the stability and 
structure of graphene but also of electronic transport. The fluctuations of normals leads to a 
modulation of the hopping integrals and are bound to affect the electronic structure [23, 121] ■ 
Knowledge of the normal-normal correlation functions is necessary for the calculation of the 
electron scattering by ripples. 

The cleavage technique that has led to the discovery of graphene has already been applied 
to other layered materials, like BN [5], so that the investigation of structural properties of 
one atom thick layers is important for a whole new class of systems. We have found that 
even fluctuations at the scale of tens of interatomic distances cannot be described by con- 
tinuum medium theory. It will be very instructive to carry out systematic experimental and 
theoretical investigations of other two-dimensional crystals to understand which properties 
are common to flexible membranes and which ones are consequences of particular features 
of the chemical bonding and interatomic interactions. 
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N 


P 


q 


L x (k) 


Ly(k) 


240 


10 


6 


24.59 


25.56 


960 


20 


12 


49.29 


51.12 


2160 


30 


18 


73.78 


76.68 


4860 


45 


27 


110.68 


115.02 


8640 


60 


36 


147.57 


153.36 


19940 


90 


54 


221.36 


230.04 



TABLE I: Details of the simulated samples. The initial, roughly squared, box is defined by 
(L x ,L y ) = (2p|ai|, <?|ai + 2a2|) where ai and a2 are the in-plane lattice vectors ai = ay^x, 
&2 = a\/3/2x+3/2y, x and y being cartesian unit vectors. 
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FIG. 1: Snapshot of the N = 8640 sample at T = 300 K. The red arrows are 70 A long. 
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FIG. 2: Sketch of a flexible membrane (solid line), h is the out of plane deviation with respect 
to the z = 0-plane (dashed line) defined by the center of mass. The unit vector ft and uq are the 
normals to each point in the membrane and in the reference plane respectively. 
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FIG. 3: Normal-normal correlation function G{q)/N for all studied samples. The solid straight 
line gives the harmonic power law behaviour Go(q)/N with k=1.1 eV. Deviations from harmonic 
behavior occur for q close to the Bragg peaks at q = 47r/3a = 2.94 A -1 and q = 47r/y(3)a = 5.11 
A -1 with a = 1.42 A and at small q where the peak of G(q) at q ~ 0.1 signals a preferred length 
scale of about 70 A. 
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FIG. 4: Normal-normal correlation function G{q)/N for the N = 8640 sample at T = 300 K and 
T = 3500 K. The solid straight line gives the harmonic power law behaviour Go{q) with k=1.1 eV 
at T = 300 and k=2.0 eV at T = 3500 K. 
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FIG. 5: Radial distribution function for the N = 8640 sample at T = 300 K and T = 3500 K as 
a function of interatomic distances in A. The arrows indicate the length of double (r = 1.31 A), 
conjugated (r = 1.42 A) and single (r = 1.54 A) bonds. 
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FIG. 6: Portion of one typical configuration of the N = 8640 sample at T = 300 K. The numbers 
indicate the bond length in A. Notice that often one of the bonds with first neighbours is much 
shorter than the other two. 
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